It is of great significance to study the bubbly medium in various fields.^[1] For example, in medicine, bubble contrast agent is used to improve the physiological tissue imaging in the field of medical ultrasound.^[2,3] In industry, bubble monitoring equipment serves as detecting the leakage in the underwater oil and gas field. In underwater acoustics, bubble film can apparently improve the conversion efficiency of the nonlinear parameter array.^[4] In ship navigation, bubble analysis can help identify the ship wake, etc.

The size, distribution, and concentration of bubbles in a gas–liquid mixture will affect not only the physical properties of the mixture, but also the nonlinear characteristics of the mixture when the bubbles are excited by sound waves.^[5] Therefore, it is necessary to obtain the size distribution of bubbles in liquid when using or avoiding the influence of bubbles.^[6] However, the bubbles usually have small sizes and most of them are in motion state, it is difficult to obtain their distribution by direct measurement. It is necessary to obtain the parameters of bubbles by other ways. Generally, the distribution of bubbles can be retrieved by acoustic and optical methods. The optical methods mostly obtain the scales through photographing or imaging, with high precision but high cost; the acoustic method has the advantages of simple measurement, strong applicability, wide measurement range, and low cost, which is more suitable for being used in complex environment.^[7,8]

In the traditional acoustic measurement methods of bubble coefficients the sound attenuation and sound velocity dispersion of bubbles are mainly used as the characteristics of sound scattering.^[9,10] The emitted sound wave will vibrate the bubbles and emits the sound energy outwards through the bubble layer. The relationship between the received acoustic data and the bubble distribution can be obtained to obtain the bubble distribution characteristics.^[11,12] The previous methods are based on the linear theory,^[1,13] so, one of their major weak points is that they are incompetent when there are scatterers other than bubbles in the ocean. Besides, it is difficult to distinguish the bubble scattering signal from the reflection or scattering signal of other objects near boundaries. Anyway, there is no effective method to filter out the influence of other factors from bubble signal.

Because the influence of other scatterers or boundaries other than bubbles corresponds to the linear scattering cross section, the nonlinear effect of bubbles makes it possible to solve the above problems to the greatest extent.^[13–17] Compared with pure water medium, bubble medium has significant nonlinear effects. These effects are independent of other scatterers or boundaries. Therefore, it is more accurate to use the nonlinear effects for the bubble distribution inversion.^[16] Most of the studies of the nonlinear scattering cross-section method do not consider the influence of non-resonant bubbles,^[18,19] or the transfer matrix with good condition number and good robustness cannot be obtained after the integral discretization, so the nonlinear inversion has not been widely used yet.^[20–22]

The contribution of this paper is that the nonlinear coefficients are selected as the basis of inversion based on the unique nonlinear acoustic effect of bubbles. Due to the good condition number of transfer matrix established by nonlinear coefficients in difference frequency by bubbles, the inverse result of size distribution of bubble groups can be obtained to be more robust. This nonlinear coefficient method seems very attractive to the ocean application because it provides high selectivity of bubbles from the other types of ocean scatterers which do not have so large nonlinear response as bubbles.

In the bubbly water medium, it is assumed that: (i) the bubble is a sphere with radius R and its vibration is spherically symmetric; (ii) the vibration of each bubble is independent (multiple scattering is not considered); (iii) the bubble size is small enough and far less than the wave length of the sound wave, so that the movement of the bubble can change at any time.

The relationship between the volume of a single bubble and the sound pressure can be described by the Rayleigh–Plesset equation

According to the relationship between the change of bubble radius and the change of volume u, the dynamic equation of the bubble volume, where the terms are kept only to the second-order terms, can be written as^[23]

Perturbative expansion of sound pressure and bubble volume is as follows:

In terms of and p₂, u₂ can be expressed as

When the bubble equilibrium radius a satisfies the distribution N(a), β(ω) can be written as

To stimulate the bubble size parameters, the nonlinear coefficient β is needed. Generally, β can be obtained indirectly by measuring B / A. The B / A is defined as the ratio of the second-order coefficient to the first-order coefficient in the Taylor expansion of the equation of state, and expressed as

The most commonly used method is the finite amplitude method, which is to measure the second harmonic amplitude due to nonlinear propagation in the medium in order to obtain B / A.

Based on the research in Ref. [8], when a finite amplitude plane wave propagates in an ideal fluid medium without viscous loss, the fundamental and second harmonic sound pressures at the angular frequency are measured respectively considering the one-dimensional sound field propagating along the x direction. The relationship between the second harmonic and the fundamental amplitude is given as follows:^[8]

To ensure the measurement accuracy, it is necessary to make the second harmonic sound pressure have not less than 6-dB signal-to-noise ratio, and determine the appropriate pulse emission intensity according to this condition.^[24] In order to avoid the influence of ultrasonic cavitation on the number of bubbles, the excessive excitation sound pressure should not be used.

Then the nonlinear coefficients are calculated according to p₁ (ω) and p₂(ω)

For example, the concerned bubble radius range is 50 μm–300 μm, the corresponding resonance frequency range is about 6.8 kHz–68 kHz, and generally the excitation frequency of 0.6 kHz–680 kHz can obtain better inversion results.

According to the measured amplitude of fundamental and second harmonic sound pressure, the measurement vector β as (M × 1)-dimensional vector of nonlinear coefficients is given by

It is assumed that the bubble radius distribution obeys the power law (PL) distribution (the power index is –4). Similarly, change the assumption that the bubble size distribution follows the Gaussian distribution (G distribution for short, variance is 1, maximum value is 5.5× 10⁹), other coefficient conditions remain unchanged, and the nonlinear coefficient is shown in Fig. 1.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Nonlinear coefficients under assumed distribution.

Figures 1(a) and 1(b) are the nonlinear coefficients, respectively, on the assumption that the distribution is G distribution and PL distribution.

The correlation equation between the bubble distribution and nonlinear coefficients is established as follows:

According to the derivation of Eq. (8), the transfer matrix K of the nonlinear parameter method (NC inversion for short) can be established as follows:

To obtain the distribution function, the transfer matrix K needs to be inversed. In the inversion problems, the condition number is often used to measure the stability of the transfer matrix. In mathematical expression, it is the product of matrix norm and matrix inverse norm. The smaller the condition number, the higher the nonsingularity of the matrix is, while the matrix is in a better state. The larger the condition number, the worse the nonsingularity of the matrix is, and the matrix is more ill-conditioned. In the numerical analysis, the influence of matrix disturbance on the eigenvalue of a given matrix is often discussed. The condition number can measure the deviation degree of matrix eigenvalue after being disturbed, and it is also an important sign to measure whether the inversion problem of transfer matrix is good or not. Therefore, the estimation of matrix condition number is of great significance for studying various matrix problems. Therefore, the condition numbers of transfer function of three nonlinear inversion methods are compared as indicated in Table 1.

Table 1.

Table 1.

Table 1. Comparison of condition number among three inversion methods. . DF inversion SH inversion NC inversion Case I 3.9602 ×1029 4.7111 × 1033 6.2547 × 1024 Case II 1.3283 × 1025 8.7555 × 1029 1.3440 × 1021

Table 1. Comparison of condition number among three inversion methods. .

In case I, the frequency range is 0.01 kHz–10000 kHz, in case II, the frequency range is 0.1 kHz–1000 kHz, and the radius range and number of points are consistent. It can be seen that the transfer matrix condition number of NC method is significantly better than that of other two methods. In the later inversion, it is more likely to show good inversion results. The is the inversion result of bubble size distribution function. With the good condition number of the transfer matrix K in the correlation equation, N and can be directly obtained by the least square estimation:

Combined with the specifically assumed distribution, the inversion analysis is carried out. The nonlinear coefficient vector ω is obtained by solving Eqs. (11) and (12). The transfer matrix K is constructed and is obtained by least square inversion. The inversion results on the assumption that the distribution obeys PL distribution and G distribution are shown in Figs. 2(a) and 2(b).

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Inversion results of nonlinear coefficient method.

It is shown that the inversion effect of nonlinear coefficient on bubble coefficient is robost without considering errors. There are some small fluctuation errors near the large radius (500 μm), which may be caused by the insufficient number of matrix points in the process of inversion, which can be improved by increasing the number of points, optimizing or interpolation. However, in order to obtain the efficiency of the inversion process, considering the accurate inversion of other radius range, such an error is acceptable.

The error analysis of the nonlinear coefficient method is carried out. The Gaussian random errors of 1%, 3%, 5%, and 10% times the standard nonlinear coefficient mean are added into the objective function of nonlinear coefficient. Similarly, 100 Monte Carlo experiments are carried out. With the B-spline interpolation, the inversion results are as shown in Fig. 3.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Error analysis of inversion results from nonlinear coefficient method: [(a)–(d)] 1%-G, 3%-G, 5%-G, and 10%-G distributions, respectively; [(e)–(h)] 1%-PL, 3%-PL, 5%-PL, and 10%-PL distributions, respectively.

Table 2.

Table 2.

Table 2. Comparison of correlation coefficient between inversion and assumed distribution after adding errors under NC inversion. . 1% 3% 5% 10% Assumed as G distribution 0.9802 0.9675 0.9420 0.9376 Assumed as PL distribution 0.8174 0.7942 0.7962 0.7794

Table 2. Comparison of correlation coefficient between inversion and assumed distribution after adding errors under NC inversion. .

Through the comparison of correlation coefficient between the inversion results of nonlinear coefficients and the assumed distribution under different errors, it can be seen that for the assumed distribution of G distribution, when the error is less than 5%, the inversion results and the assumed distribution can have a correlation degree of more than 0.9; under the conditions of errors 5% and 10%, although the inversion results are near the assumed distribution, they fluctuate greatly at small radius, the correlation coefficient still remains high. When the assumed distribution obeys the PL distribution, for the number of bubbles changing greatly in the small radius range, the inversion effect is relatively unsatisfying, which is reflected in the case of larger difference between the blue line and red line in the small radius and the smaller correlation coefficient; in most of radius range, the inversion effect is ideal.

The influence and effect of different distributions. It can be seen that the effect of NC inversion is different when the assumed distribution is different. If the distribution is assumed to follow the G distribution, the number of bubbles changes less with the bubble radius, so the inversion effect is better; if the distribution is assumed to follow the PL distribution, the number of bubbles changes more in the small radius range, so the inversion does not achieve good results. The correlation coefficient of the inversion curve and the assumed distribution curve are smaller than that of the G distribution, and the effect is not so stable as the G distribution. In the range of more concentrated bubbles, in order to have a good inversion effect when the change is severe, we can optimize further the interpolation algorithm to improve. There is no further discussion here. Using a similar matrix equation inversion method to invert the bubble distribution, the correlation coefficients between the assumed distribution and the inversion results with different errors are obtained.

The comparioson between the results from other two nonlinear inversion methods are shown in Table 3. The specific inversion results are shown in Fig. 4. Apparently, based on the same optimization method, NC inversion shows better inversion effect and robustness compared with other two inversion methods.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Error analyses of inversion results of DF inversion method: [(a)–(d)] 1%-G, 3%-G, 5%-G, and 10%-G distributions, respectively; [(e)–(h)] 1%-PL, 3%-PL, 5%-PL, and 10%-PL distributions, respectively.

Table 3.

Table 3.

Table 3. Comparison of correlation coefficient between inversion and assumed distribution after adding errors under DF and SH inversion. . Inversion methods Errors\Assumed 1% 3% 5% 10% DF inversion G distribution 0.6153 0.5720 0.5997 0.6077 PL distribution 0.7189 0.7744 0.7917 0.6093 SH inversion G distribution 0.3245 0.2021 0.2081 0.3530 PL distribution 0.1048 0.0355 0.0642 0.0892

Table 3. Comparison of correlation coefficient between inversion and assumed distribution after adding errors under DF and SH inversion. .

To test the NC inversion’s practicability and effectiveness, laboratory experiment is carried out. The process of inversion of bubble distribution by nonlinear coefficient method is simulated through specific measurement as follows.

The projecting transducer and the receiving hydrophone are respectively placed on both sides of the inhomogeneous mixed medium sample with thickness. The projecting transducer emits single frequency pulse to the bubble water medium to be tested at x = 0. The receiving hydrophone receives the pulsed sound pressure signal along the axial direction. The far-field conditions of the transceiver are satisfied as shown in Fig. 6.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Error analyses of inversion results of SH inversion method: [(a)–(d)] 1%-G, 3%-G, 5%-G, and 10%-G distributions, respectively; [(e)–(h)] 1%-PL, 3%-PL, 5%-PL, and 10%-PL distributions, respectively.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Experimental settings of measuring nonlinear coefficients of bubble medium.

In the experimental water tank, the combined projecting transducer is placed on one side of the tank at the position of x = 0, and the receiving hydrophone is placed on the other side of the tank at the position of x = R, where settings meet the far-field condition. The combined projecting transducer and the hydrophone are set to be coaxial. In Fig. 6, 1 is connected to the transmitting equipment; 2 is connected to the receiving equipment; 3 represents the combined projecting transducer; 4 represents the receiving hydrophone; 5 is a bubble electrolyzer; 6 is the produced bubble group; 7 is the water tank; 8 is the pure water medium; 9 is the coordinate axis, representing the acoustic axis.

Because the measurement system needs to cover a large frequency range, the multiple transmitter transducer combinations are actually used in the measurement. The center frequencies are 90 kHz, 150 kHz, 250 kHz, 350 kHz, 500 kHz, and 1000 kHz respectively. All the transmitting transducers are of circular planar piston type, with diameters in range of 70 mm–200 mm. The size of the water tank used in the experiment is 1200 mm × 780 mm× 1000 mm. The center frequency of the hydrophone is 250 kHz, the sensitivity is –214.3 dB, the size of the probe is about 20 mm, and the receiving frequency range is 10 kHz–5 MHz, which can maintain a relatively stable electrical signal output. The experimental device is symmetrically placed in the water tank in terms of depth and width, the bubble electrolyzer is placed in the center of the water tank, and the transceiver is symmetrically placed on both sides with a spacing of 1 m. For the frequency range of transmission, the layout meets the far-field conditions, and does not consider the case that the transceiver signal is non-plane wave. Meanwhile, the sizes of the transducers and hydrophone can also be ignored. The experimental system block diagram is shown in Fig. 7.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Experimental system block diagram.

Coefficients are set be as follows: water medium density ρ₀ = 1000 kg / m³, and sound speed c₀ = 1500 m/s. Based on the prior knowledge, considering the radii of bubbles produced by electrolysis are in a range of 15 μm–120 μm, the excitation frequency range is chosen as 60 kHz–1200 kHz. Assuming that the number of radius points N = 1000, the number of frequency points M = 1000, damping coefficient δ = 0.1, specific heat ratio γ = 1.4, and hydrostatic pressure P₀ = 1.088× 10⁵ Pa. In order to reduce the random error caused by the instability of bubbles in the experiment, the number of measurement pulses is 10, and the relatively stable measurement results are obtained by taking the average value.

By measuring the sound pressure and interval between the fundamental and the second harmonic of the received signal, the values of nonlinear coefficient β of the medium at different excitation frequencies can be obtained from Eq. (12).

According to the measured frequency range and frequency sampling interval, a nonlinear transfer matrix is established, and the measured nonlinear coefficients are sorted into an objective function matrix. The inversion results optimized by B-spline interpolation and regularization are as follows.

According to the same measurement data, we can also use the SH inversion method to invert and obtain the following results.

Compared Fig. 9 with Fig. 10, NC inversion method shows a reasonable bubble distribution, while SH inversion method results in an obviously divergent result under the same fittings and optimization methods (The tail peak of the line is warped, which does not conform to the general law of bubble distribution, so the inversion results are divergent). It reflects the better robustness of NC inversion method. As a new bubble distribution inversion method, the NC inversion may be widely used in acoustical measurement in the future.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Nonlinear coefficients measured in experiment.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Results of bubble size distribution retrieved by NC inversion method.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. Results of bubble size distribution retrieved by SH inversion method.

According to the frequency range and measurement points selected in the experiment, the inversion effect of NC inversion method can be estimated. When 5% error is introduced, the inversion effect is shown in Fig. 11, with the bubble group assumend to have the G distribution.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. Expected inversion effect under experimental conditions.

It can be seen that under the same assumption error, the inversion effect is not ideal for the frequency range under experimental conditions. Next, the inversion effects in different frequency ranges are discussed below.

The nonlinear coefficient data of six frequency ranges are selected for inversion, and the results are shown in Fig. 12. The selected frequency ranges are: 60 kHz–400 kHz, 60 kHz–600 kHz, 60 kHz–800 kHz, 200 kHz–600 kHz, 200 kHz–1000 kHz, and 600 kHz–1000 kHz. Respectively corresponding to Figs. 12(a)–12(f).

Fig. 12.

	Figure Option ViewDownloadNew Window
	Fig. 12. Inversion effects of β data in different frequency ranges.

The frequency range discussed is based on the resonance frequency of the concerning bubble radius and the limitations of the transceiver. The frequencies 1000 kHz and 60 kHz are the upper and lower limits of emission frequency, while 400 kHz is the upper limit of the resonance frequency concerning bubble radius. The frequencies 200 kHz, 600 kHz, and 800 kHz are 0.5, 1.5, and 2 times the upper limit of the resonance frequency respectively.

The loss of nonlinear coefficient in the low frequency band will result in the confusion of the whole order of magnitude. The loss of high frequency nonlinear coefficient will lead the curve to sharply change, which is reflected in the displacement of the concave position to the small radius. From this point of view, the objective function information in low frequency plays a more important role in realizing the inversion. Through the study of the objective function in the theoretical derivation stage, it can be seen that the nonlinear coefficient has more fluctuations in the low-frequency segment, which means that the information in the low-frequency segment has a greater influence on the inversion process. Therefore, when using NC inversion method, the abundant low-frequency measurement results are needed for the accuracy of inversion.

The nonlinear coefficient inversion method takes the nonlinear coefficient as the sound field information to retrieve the bubble group size distribution and directly obtain the physical coefficients of the bubble without using the acoustic coefficients of the bubble. The main advantages are indicated below.

(I) For nonlinear inversion methods to obtain the coefficients of the bubble, the inversion error caused by the influence of boundary or other scatterers on the acoustic field information can be stripped.

(II) In the nonlinear inversion methods mentioned above, the matrix condition number obtained by kernel function discretization is not robust enough. The inversion result can be obtained by simple least square method. The calculation is simple and effective, and the inversion result is better. But their robustness degrees are different from each other.

(III) The robustness degree of NC inversion is better than those of other two nonlinear inversion methods. This conclusion is also verified in experiment. The NC inversion reaches an more stable result than SH inversion.

Meanwhile, when the error is introduced, the nonlinear coefficient inversion method does not show better robustness in some bubble ranges nor large error, which leaves a further research of nonlinear inversion methods, especially the nonlinear coefficient inversion method.

(IV) In the experimental measurement, the measurement of nonlinear coefficient in the lower frequency band has an influence on the order of magnitude of inversion, and the measurement of nonlinear coefficient in the higher frequency band has an influence on the shape of inversion curve. Therefore, it is more important to obtain the abundant low-frequency measurement data.